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We develop an explicit description of a time-dependent response of fermionic condensates to 
perturbations. The dynamics of Cooper pairs at times shorter than the energy relaxation time can 
be described by the BCS model. We obtain a general explicit solution for the dynamics of the BCS 
model. We also solve a closely related dynamical problem - the central spin model, which describes a 
localized spin coupled to a "spin bath" . Here, we focus on presenting the solution and describing its 
general properties, but also mention some applications, e.g. to nonstationary pairing in cold Fermi 
gases and to the issue of electron spin decoherence in quantum dots. A typical dynamics of the 
BCS and central spin models is quasi-periodic with a large number of frequencies and stable under 
small perturbations. We show that for certain special initial conditions the number of frequencies 
decreases and the solution simplifies. In particular, periodic solutions correspond to the ground 
state and excitations of the BCS model. 



I. INTRODUCTION 

Recent experimentsi on cold fermion pairing in the 
vicinity of a Feshbach resonance offer a unique opportu- 
nity to control the strength of pairing interactions be- 
tween fermions by a magnetic field. This opens up an 
exciting possibility to explore fundamentally new aspects 
of the nonequilibrium pairing following an abrupt change 
of the coupling strength. 

Interestingly, the dynamical pairing in this regime can 
be linked to a seemingly unrelated spin model. The 
latter describes the interaction of a localized (central) 
spin with a "spin bath" of environmental spins. The 
central spin model has reemerged recently in the con- 
text of experiments^ on electron spin dynamics due to 
the hyperfine interaction with nuclei in GaAs quantum 
dots. Since single electron quantum dots are now experi- 
mentally accessible^, and are considered one of the most 
promising candidates for solid state qubits, it is also im- 
portant to understand the effect of this interaction on the 
coherence of the electron spin. 

In this paper, we present a theory that provides an 
accurate description for a range of phenomena in the dy- 
namical pairing and central spin problems. As we will see 
below, there is an intimate connection between these two 
problems that will enable us to treat them on an equal 
footing. 

The study of the dynamics of the superconducting 
state in metals founded on microscopic principles has 
begun during the decade following the advenlii of the 
BCS theory (see Ref. i for a review). The simplest 
theory for nonstationary processes in superconductors is 
based on the time-dependent Ginzburg-Landau (TDGL) 
equation-'^-- for the order parameter A(t). The TDGL 
approach is valid when quasiparticles are able to reach a 
local equilibrium quickly on the characteristic time scale 
of the order parameter variation ta (— at zero tem- 



perature) . This requirement usually limits the applicabil- 
ity of the TDGL theory to situations where mechanisms 
destroying Cooper pairs are effective, such as a narrow 
vicinity of Tc or a large concentration of magnetic impu- 
rities. 

An alternative theorjiSiiS is the Boltzmann kinetic ap- 
proach, which describes the dynamics in terms of a ki- 
netic equation for the quasiparticle distribution function 
coupled to a self-consistent equation for A(i). However, 
this scheme is justified only when external parameters 
change slowly on the ta time scale. 

As pointed out in Ref. 11, cold fermionic gases can be 
in a regime where the notion of the excitation spectrum is 
irrelevant and neither TDGL nor the Boltzmann kinetic 
equations are valid. Indeed, in these systems external 
parameters such as the position of a Feshbach resonance 
can change on a time scale tq <C ta, Tc, where is the 
quasiparticle energy relaxation time. On the other hand, 
the energy relaxation is slow, while the lifetime of the 
samples is limited. It is therefore desirable to develop 
a theory that describes nonstationary pairing effects in 
this regime. 

The nonequilibrium Cooper pairing at times t ^ is 
non-dissipative and in a translationally invariant system 
can be described by the reduced BCS model 
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It is interesting to note that the use of this model for de- 
scribing the dynamics of a homogeneous superconducting 
state in metals in the non-dissipative regime, t <^ t^, has 
been justified^ quite generally starting from nonstation- 
ary Eliashberg equations. 

Here, we are interested in a situation when a system 
at zero temperature is out of an equilibrium at t = 0. 
The goal is to determine the subsequent evolution of the 
initial state. In particular, this includes the case when at 
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t = the coupling constant has been abruptly changed 
from g' to g. 

In the absence of translational invariance, e.g. in a 
dirty superconductor or in a finite system, Hamiltonian 
has to be appropriately generalized^^. First, let us 
discuss the familiar case of electron-electron interactions. 
If there are no spatial symmetries, but the time-reversal 
invariance is preserved, single-particle orbitals ej are de- 
generate only with respect to spin. For each orbital ej 
there is a pair of states, \j t) and \j |), related by time 
reversal symmetry. The pairing occurs between time re- 
versed states, i.e. 

Hbcs = ^Aoh'^ " ^T^i'^'JiCgT (1-2) 

Deviations of the coupling strength g from a constant in 
q and j are small and sample-dependent 

For a two species Fermi system, such as a two- 
component mixture of fermionic atoms in a trap, spin 
up and down states in (|1.2I) have to be identified with 
the two species of fermions. With this replacement, 
argumentsi^iii^^ leading to Hamiltonian H1.2|l are valid 
in the weak coupling regime. 

Hamiltonian is obtained from (|1.2|) by specializing 
to a translationally invariant situation. In this case, time 
reversed states are |p t) and | — p i). One can also take 
a thermodynamic limit by taking the number of orbital 
levels to infinity. However, in view of applications where 
the effective number of levels is finite and even small (see 
below), we deal with the general Hamiltonian (|1.2|l here. 

Remarkably, the BCS model turns out to be closely 
related to a model describing the Heisenberg exchange 
interaction of a single localized spin (the central spin) 
with a number of environmental spins. The Hamiltonian 
reads 

n-l 

ffo = 51 • ~ (1-3) 

where Kg is the central spin, Kj are environmental spins, 
7j are (nonuniform) coupling constants, and B represents 
the magnetic field^^. 

The central spin model has an interesting history of its 
own. For example, it emerged in the studies of electron 
spin dynamics in disordered insulatorsii and of coher- 
ent spin tunnelling in ferromagnetic grainsi^. More re- 
cently, it attracted considerable attention as a model for 
decoherence of qubitsi^iSfliSiiS^. In particular, it can be 
used to model^^-^° the hyperfine interaction of a localized 
electron spin with nuclear spins and the spin-dependent 
transport in GaAs quantum dots. In these cases, it 
provides an adequate description of spin dynamics at 
times shorter than the nuclear spin relaxation time in 
the regime where the orbital level spacing in the dot is 
much larger than the temperature and typical interaction 
energies. 

The main result of this paper is an explicit general so- 
lution for the dynamics of the BCS and central spin mod- 



els. In particular, we determine as functions of time the 
order parameter A(t) and the expectation value of the 
central spin (equations l|3.22|l and (|3.24|) ') as well as the 
dynamics of the remaining degrees of freedom for arbi- 
trary initial conditions. Here, we concentrate on present- 
ing the solution and describing its general properties. We 
also mention several applications to specific issues such 
as nonequilibrium pairing and that of decoherence due to 
the hyperfine interaction, but leave a detailed discussion 
for the future. 

The solution is based on the integrability2^*SS*2& of the 
BCS and central spin models. This important property 
has been largely underestimated in part because it was 
discovered outside the main physical context of these 
models. Besides, due to the infinite range of interactions 
in Hamiltonians H1.2I) and 1)1. the mean field approxi- 
mation is exac1iS2iS£ for these models in the limit of large 
number of particles. In this approximation, the BCS and 
central spin models can be mapped onto a classical non- 
linear system (see below). However, while the mean field 
simplifies the description of equilibrium properties to an 
extent where no advanced techniques are required, it is 
the integrability of the resulting classical system that en- 
ables us to solve dynamical problems. The BCS solution^ 
for the ground state and excitations of the BCS model is 
recovered as a periodic case of the general solution (see 
the discussion following (|2.6|) and Section llV|) . It is also 
interesting to note that, as we will see, models (|1.2|l and 
(|1.3|) . being in many respects classical, have distinct ro- 
bust features that are preserved when the integrability is 
destroyed. 

The Dicke model'^^ that describes an ensemble of two- 
level systems coupled to a bosonic mode is also within the 
scope of our construction. In fact, it belongs to the same 
class^^ of integrable models as the BCS and central spin 
models. The Dicke model has been recently adopted— 
to the study of a dynamical coexistence of atoms and 
molecules in cold Fermi gases. In the strong coupling 
regime, it seems to be an appropriate alternative to the 
BCS model for these systems. 

The general nonstationary solution of the BCS and 
central spin models is in terms of hyperelliptic func- 
tions - multiple variable generalizations of ordinary el- 
liptic functions. These functions frequently arise as so- 
lutions of integrable equations and have well known an- 
alytical properties'^. 

We show that for most initial conditions the dynamics 
is gwasi-periodic with a number of independent frequen- 
cies equal to the number of degrees of freedom, n, and 
evaluate the frequencies in terms of the integrals of mo- 
tion. The typical motion uniformly explores an invariant 
torus - an n-dimensional subspace of the 2n-dimensional 
phase space allowed by the conservation laws. These fea- 
tures of the typical motion are stable against small per- 
turbations that destroy integrability. 

Further, we identify certain special values of integrals 
of motion for which the number of independent frequen- 
cies, m, becomes less than the number of degrees of free- 
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dom. For these degenerate cases, we were able to explic- 
itly reduce the motion of the BCS and central spin models 
with n degrees of freedom to that of same systems with 
only m degrees of freedom. The solution progressively 
simplifies as the number of independent frequencies de- 
creases. For example, solutions characterized by a single 
frequency (m = 1), i.e. periodic trajectories, are given 
by trigonometric functions; degenerate solutions with two 
independent frequencies (m = 2) are given by a combi- 
natio n of t rigonometric and elliptic functions (see also 
Refs. Illisol and the end of Section Wv} etc. 

Periodic solutions {m ~ I) occupy a special place 
among degenerate solutions. As we show in Section IWI 
they reproduce the BCS solution for the ground state and 
excitations of the BCS model. 

However, it should be emphasized that, unlike the gen- 
eral solution, degenerate solutions are nonrepresentative 
of the dynamics of the BCS and central spin models and, 
with the exception of the ground state, are expected to 
be unstable. 

The paper is organized as follows. In Section|nl we use 
a mean field approximation to map the BCS and central 
spin models onto equivalent classical models. Section UTTI 
contains an explicit general solution for the dynamics of 
both problems. In Section IIVI we consider degenerate 
solutions. Section discusses some open problems and 
possible applications of our results. 



II. MEAN FIELD APPROXIMATION 



Our goal is to determine the time evolution accord- 
ing to Hamiltonian (|2.1|l of an arbitrary initial distribu- 
tion of pseudospins. The mean field approximation con- 
sists in the replacement of the effective field seen by each 
pseudospin in BCS Hamiltonian (12.11) with its quantum 
mechanical average, hj{t) = {—2Ax{t),—2Ay{t),2ej), 
where A{t) is the BCS gap function 

A(t) ^ A,{t) - ^Ay{t) ^ 9j2{K-{t)) (2.3) 



In this approximation, each spin evolves in the self- 
consistent field created by other spins 



BCS, 



(2.4) 



Since equations 1)2. 4|l are linear in Kj, we can take their 
quantum mechanical average with respect to the time- 
dependent state of the system to obtain 



Sj hjxsj 



^3 ^ i-^aJx^-gJy^'^.ej) J^^s, (2.5) 



where Sj{t) — 2{'Kj(t)). Evolution equations H2.5|l con- 
serve the square of the average for each spin, i.e. = 



const. If spins initially were in a product statc'^^, — 1. 

Note from equation (|2.2|l the following correspondence 
between components of Sj (t) and the normal and anoma- 
lous (Keldysh) Green functions at coinciding times: 



Our starting point is the mean field approximation, 
which enables the mapping of the BCS and central spin 
models onto equivalent classical nonlinear systems. We 
also derive the equations of motion and describe the re- 
lationship between the BCS and central spin models. 

The discussion of the mean field is facilitated by rep- 
resenting BCS Hamiltonian (|1.2|l in terms of Anderson 
pseudospin- 1/2 operators2I. 



H 



n-l 



n-1 



BCS = 2^ 2eji^J - g L+L. 

j=0 g=0 



L = ^K, (2.1) 



where n is the number of single-particle orbitals. Pseu- 
dospin operators are related to fermion creation and an- 
nihilation operators via 



(2.2) 



Fj{t) = -i{[c3^it),c,^it)]) ^ isjit) 

Equations H2.5|l were derivedsi? for phonon supercon- 
ductors within the general framework of nonstationary 
Eliashberg theory in the coUisionless regime t <C t^. A 
linearized version of these equations was considered in 
Refs and m 

Due to the infinite range of interactions between spins 
in H2.1(l . the mean field approximation is exact in the 
thermodynamic limit. For a system with a finite number, 
N , of particles (spins) we expect, based on an analysis of 
leading finite size corrections^ to the mean field, equa- 
tions (|2.5I) to be accurate for large A'^ at times t < t*N^, 
where t* and 77 > do not depend on N. 

We see that Sj — 2{Kj) have all properties of classical 
spins governed by a Hamiltonian 



BCS 



(2.6) 



Pseudospins are defined on unoccupied and doubly occu- 
pied pairs of states \j t) and \j |), where they have all 
properties of spin-1/2. Singly occupied pairs of states, 
on the other hand, do not participate in pair scattering 
and are decoupled from the dynamics. 



and usual angular momentum Poisson brackets. Thus, 
the problem reduces to determining the time evolution 
of n classical spins according to Hamiltonian (|2.6|l . 

The BCS solution for the ground state corresponds to 
the minimum of H2.6|l and is obtained by aligning each 
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spin in H2.6|l antiparallel to the field acting upon it. A 
pair excitation'^'' (excitation of the condensate) of energy 
(ej — + Aq, where pi and Aq are the chemical po- 
tential and the equilibrium gap respectively, is obtained 
by flipping the spin Sj. The C/(l) symmetry of the BCS 
order parameter A is equivalent to the symmetry of (|2.t)|) 
with respect to uniform rotations of all spins around the 
z-axis. Due to this symmetry, spin configurations that 
determine the ground state and excitations can rotate 
around the 2;-axis with a frequency 2/i at no energy cost. 
In the presence of a particle-hole symmetry /i = and 
these configurations are stationary with respect to Hamil- 
tonian (|2.t)|l (see Section IIVI where we recover the BCS 
solution as a periodic case of the general time-dependent 
solution). 

Next, we turn to the central spin model and its rela- 
tionship to the BCS model. The central spin Hamiltonian 
turns out to be a member of an integrable family of 
Gaudin magnet^^ - n Hamiltonians of the form: 

il^ = 2Y^ ' - akl 9 = 0,. ..,71-1 (2.7) 

3=0 ^« ~ 

where Cj and a are arbitrary parameters. All Hamil- 
tonians Hq commute with each other and with the z- 
component of the total spin, Hq, for arbitrary 

tj and a. The central spin model H1.3|l coincides with _ffo 
in equation (|2.7|l if we choose eg — tj = "ij^j and a = B. 

Note that if a number of orbitals in 1)2.1(1 are degen- 
erate, the magnitude of their total spin X^e =const 
conserved by Hamiltonian H2.1|l . In this case one can re- 
place Kj — > X^e =const -^j Hamiltonians H2.1|l and (|2.7|l 
and sum over nondegenerate orbitals only. 

Now consider the following linear combination of Hq. 



Hamiltonian H1.3(l is 



q=0 



.9=0 



2 

a 



const 



We see that for a — 2/g the expression in the square 
brackets coincides with BCS Hamiltonian (|2.1|l . There- 
fore, for a = 2/5, the BCS Hamiltonian commutes with 
all Hq and thus belongs to the same class of integrable 
modelsS^. 

In the same way as we did for the BCS model, 
we employ the mean field approximation to derive 
classical Hamiltonians that govern the evolution of 
Sj{t) = 2(Kj(t)) for Hamiltonians lfT7|) : 



Hq = Y.'^^ 

3=0 ^« 

hbcs = J2 
3=0 



as„ 



J+J- 



(2.8) 



3=0 



(2.9) 



The validity of the mean field approximation for central 
spin model 1)1.3(1 is subject to similar considerations as 
for the BCS model (see the second paragraph following 
equation 1(2.5(1 and also Ref . JOJ . 

Thus, the problem of determining the dynamics of the 
BCS and central spin models reduces to solving equations 
of motion for Hamiltonians Hq and Hbcs in system (12.811 . 
To obtain solutions for BCS 1(1.2(1 and central spin 11. 3() 
models, one has to choose 



a = - BCS 
9 



(2.10) 



eo - 



= 7j a = B central spin 



In particular, the classical counterpart of central spin 



Hamiltonians 1(2.8(1 Poisson-commute with each other 
and with the 2;-component of the total spin cx Hq. 
Each Hamiltonian in system ((2.8(1 describes an evolution 
of n spins in an 2n-dimensional phase space (two an- 
gles for each spin) and has n integrals of motion (the en- 
ergy and the remaining Hamiltonians from system ((2.8(1 ) . 
Therefore, all Hamiltonians 1(2.8(1 are classical integrable 
models. 



III. THE SOLUTION 

In the previous section we mapped the BCS and cen- 
tral spin models onto a classical integrable system ((2.8(1 . 
Nevertheless, even though Liouville's theorem^Si guaran- 
tees a formal integrability in quadratures of classical in- 
tegrable models, an explicit solution for the evolution of 
the original dynamical variables is not always possible. 
Fortunately, this turns out to be not the case for the 
BCS and central spin problems. For these models, as de- 
tailed below, we were able to obtain an explicit general 
solution of equations of motion. 

Before we proceed, let us discuss generic featurea2& ex- 
pected of the dynamics of a classical integrable model 
with n degrees of freedom. The motion is confined by 
conservation laws to an n-dimensional subspace of the 
2ri-dimensional phase space. This subspace (invariant 
torus) is determined by initial values of integrals of mo- 
tion and is topologically equivalent to an n-dimensional 
torus. The motion on the invariant torus is characterized 
by n angle variables, (j)k and the corresponding angular 
frequencies, (pkit) — ^kt + ^fc(O). The frequencies de- 
pend only on integrals of motion, while constants (j)k (0) 
are determined by initial values of remaining degrees of 
freedom. Since for most initial conditions all frequen- 
cies ilk are independent"^^, typical trajectories uniformly 
explore the entire torus. All these properties are not 
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affected by small perturbations destroying integrability 
(Kolmogorov- Arnold-Moser theorem) . 

Now we turn to the solution for the dynamics of the 
BCS and central spin models. The solution consists of 
two main steps. The first one is a change of variables 
that casts equations of motion into the form of a known 
mathematical problem. In the second step, we use the 
solution of this problem to obtain dynamical variables 
Sj{t) = 2{'Kj{t)) as explicit functions of time. 

The change of variables is facilitated by defining a 2 x 2 
matri3*2& that depends on dynamical variables s_,- and also 
on an auxiliary parameter u 



A{u) B{u) 
" ' B*{u) -A{u) 



where matrix elements of C are 

n— 1 n— 1 



J=0 



(3.1) 



(3.2) 



The eigenvalues of this matrix, ±w(m), depend on dy- 
namical variables only through integrals of motion - Hj 



and 



2Hi 



3=0 



(3.3) 

where we defined two polynomials, Q2n{u) and Pn{u) 
that will be frequently used in subsequent calculations 

2ri-l 2n 

Q2n{u) = \{{u-E,) = Y, aku"" a2n = 1 (3.4) 

j=0 fe=0 



(3.5) 



Note that because C is Hermitian for real u, its eigen- 
values are real and therefore the polynomial Q2n{u) is 
positive definite on the real axis. We will often refer to 
the polynomial Q2n(u) as the spectral polynomial. 

Following an algebraic version of variable separation 
mcthodt^2ii£ we introduce n — 1 variables, Uj, as zeroes 
of B{u) - one of the off-diagonal matrix elements of C. 
Variables Vj, canonically conjugate to Uj, are given by 
one of the eigenvalues of matrix C at u = Uj. 



B{u,) = 



-A{u,) j = 



1 



(3.6) 



Although we will not use the Hamilton-Jacobi method, 
we remark that variables {uj , Vj ) separate^ Hamilton- 
Jacobi equations for Hamiltonians H2.8|l . Since Uj are 
zeroes of B{u), we can rewrite the matrix element B{u) 
as 



B{u) = J_ 



1X1=1 (u-Uj) 



X 



Rn-l{u) 
Pn{u) 



(3.7) 



This equation is useful for expressing original dynamical 
variables - components of spins Sj - through separation 
variables Uj and J_. For example, using equations 1)3. 2|l 
and H3.7|) . we have 



res B{u) — J. 



■^n-i(ej) 



(3.8) 



where the prime denotes a derivative with respect to u. 

Equations of motion for central spin Hamiltonian H2.9() 
in terms of new variables Uj and J_ can be derived 
by evaluating Poisson brackets between the Hamiltonian 
and matrix elements of C. We have 



iBy/Q2niUj) 



B 



(3.9) 



j_ = iBs^ = iA J_i?„_i (0) A = — [] 7j 



Similarly, one obtains equations of motion for BCS 
Hamiltonian H2.6I) 



J_ 



2^VQ2«(^3^) 

U„,^jiUj ~ Um) 

' n-1 n-1 

J=0 3 = 1 



(3.10) 



In equations H3.9|l and H3.10|l as well as in the rest of the 
paper, with no loss of generality, we shifted parameters ej 
in the BCS and central spin Hamiltonians by a constant 
so that eo = 0. 

In the reminder of this section we obtain an ex- 
plicit general solution for the expectation value of each 
(pseudo)spin, Sj{t) — 2(Kj(t)), as a function of time. 
In particular, this includes the expectation value of 
the central spin, So(i), and the BCS gap function, 
A(i) =5J_(t). 

The solution is based on the observation that, with a 
help of elementary algebra, equations of motion (|3.9|l and 
(|3.1UI) can be cast into a form of a known mathematical 
problem. Specifically, one can rewrite equations H3.9|l and 
(|3.1U|) for variables uj as 



n-l i-lj 

Uj auj 



E 



I VQ2jUjj 



— dxi I — 1, 



(3.11) 



where the polynomial Q2n{u) is defined in equation (|3.3|l . 
To obtain equations of motion (|3.1U|) for the BCS Hamil- 
tonian one has to choose 

= j(ci, . . . , c„_2, 2t + c„_i) (3.12) 

while for the central spin model one has 

= -iX{t + ci,.. . ,c„_2,c„_i) (3.13) 
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where ci, . . . , c„_i are arbitrary constants. 

Differential equations H3.11|l constitute a well-known 
mathematical problem called Jacobi's inversion problem 
(see, for example, Ref. 31 and references therein). The 
solution of equations (|3.11l) can be expressed^ through 
hyperelliptic Abelian functions - multiple variable gener- 
alizations of ordinary elliptic functions. These functions 
are often encountered as solutions of integrable equa- 
tions and have well known analytical properties. They 
are also implemented in standard mathematical software 
packages^. 

Riemann theta function of genus g (in our case 
g = n — 1) is defined as the following sum over all g- 
dimensional integer vectors 



(x|r) 



E 

mGZS 



exp [i7r(m"'"Tm -H2x^m)] (3.14) 



where r = uj'lu^^, and oj and uj' are g x g matrices to be 
specified below. Klenian a- and C-functions of genus g 
are defined through the Riemann theta function aa^ 



cr(x) Cexp[x'^77(2tj)-ix]6' ((2w)-1x|t) 



0(x) 



91n(T(x) 

dxi 



(3.15) 



In our case, the g x g matrices oj and rj (matrices of pe- 
riods) that appear in the definition of hyperelliptic func- 
tions H3.14|l and H3.15|l are 



bk \/Q2n{u) 

du 



'kl 



2n~l 

E 



(3.16) 



(fc - l)ak+iu 



k-l 



where contours of integration bk go around branch cuts 
of a hyperelliptic curve z{u) = \/ Q2n{u) as shown on 
Fig. 1. Matrices u)' and 77' are also defined by equations 
(|3.16|) with the replacement of the contours of integration 
hk hk (see Fig. 1). 

The solution^ of differential equations (|3.1HI for ar- 
bitrary initial conditions states that functions itp(x) are 
zeroes of the following polynomial of a degree g = n — 1: 



i?„_l(7/,i) -7/"-l-^/fc(x) 



fc=l 



where coefficients of Rn-i{u,t) are 

/fc(x) -a(x+d)-a(x-d) 



(3.17) 



(3.18) 



Here vector x is defined by equations (|3.12() and (|3.13() 
for the BCS and central spin problems respectively. The 
argument of ^-functions in equations H3.18|l is shifted by a 
vector of constants d that has the following components: 



Eo \/Q2n(w) 



1 = 1,. 



,71—1 




FIG. 1: Riemann surface of a hyperelliptic curve z{u) = 
-\/ Q2n{u) of genus'** g = n — 1 showing contours of inte- 
gration bk and hk with k = 1, . . . ,g. These contours appear 
in the definition of matrices of periods 13.161 of hyperellip- 
tic functions. Contours bk go around g = n — 1 branch cuts 
connecting first 2g = 2n — 2 roots Ek of the spectral poly- 
nomial Q2n{u). Integration along the contours is clockwise. 
Note that since the spectral polynomial is positively defined 
the branch cuts are parallel to the imaginary axis. 



where Eq is one of the roots of the spectral polynomial 
Q2n{u) (see Fig. I and equation (|3.4() '). Note that, ac- 
cording to its definition H3.3|) . the spectral polynomial 
Q2n(u) depends only on integrals of motion and on pa- 
rameters Ej and a. 

Finally, constants aj, in equations H3.18|l are obtained 
from the expansion 



1 



ai 



k=0 



Z 



0(1) (3.20) 



where ak are coefficients of the spectral polynomial 
Q2n{u) defined by equation (|3.4|l . For example, a„_i = 
a2n~i/2, a„_2 = a2n-2/2 - ai„_i/8 etc. 

To obtain quantum mechanical expectation values of 
original spins Sj{t) = 2(Kj(i)) for central spin model 
(|1.3|) and pseudospins (|2.2I) for BCS model H1.2|l . we ob- 
serve that, according to equation H3.8|l . they can be ex- 
pressed through the polynomial i?„_i(u,t) given by ex- 
pression l|3.17|l . We have 



-{t) = J-{t) 



(3.21) 



Recall that coefficients of the polynomial P„ (m) , defined 
by equation (|3.5I) . depend only on parameters ej. 

Similarly, one can use an expression for the matrix el- 
ement A{u) defined in equation H3.2|l in terms of vari- 
ables Uj to obtain z-components of (pseudo)spins as func- 
tions of time. However, these components are not in- 
dependent and can be also obtained from the equation 



(3 19) I j ~ where the sign of s| is determined 



by initial conditions 
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Integrating equations (|3.9() and H3.10|l for the evolution 
of X and y components of the total (pseudo)spin with the 
help of equation H3.18|l , we obtain 

J-ft)=c„e"^^' ^i'' + ^i A(i)=.9J_(i) (3.22) 
cr(x — dj 

where vector x is defined by equations (|3.12() and H3.13|l . 

the frequency (3 is different for the BCS {(3bcs) ^^^d cen- 
tral spin (/3o) models 

n-l 

/?BCS = 5-^2 + 2 ^ - a2„-i /3o = aiA (3.23) 

3=0 

and A is defined in equation H3.9|l . 

Average components of the electron (central) spin, 
so(t) = 2(Ko(i)) are obtained from equation H3.21|l (re- 
call that £0 = 0) 

s,{t) = [Ci(x - d) - Ci(x + d) - a,] ^|^±^ 

(3.24) 

where = {\cn)/B. 

To summarize, quantum mechanical averages of 
(pseudo)spins for the BCS and central spin models, in- 
cluding the electron (central) spin and the BCS gap func- 
tion A(i) are hyperelliptic functions (equations (|3.14(l 
and (|3.15|) ) of the time variable and n — 1 constants Ck 
(equations H3.21|l . H3.22|l . and (|3.24|l ). These functions 
are specified by their matrices of periods uj, uj' ^ and rj. 
The latter depend only on integrals of motion according 
to equations (I3.3|l and H3.16|l . The remaining n out of 2n 
initial conditions fix constants Ck in equations (|3.12l) and 

(Tm . 

An important characteristic of a quasi-periodic motion 
are its (Fourier) frequencies. To determine the frequen- 
cies for the BCS and central spin models, we observe from 
equation l|3.15|l that the time variable enters the solution 
only through combinations of 6*'^'=*, where 

Vlk = iT^i^k^ ^0 = Pbcs BCS 

(3.25) 

Qk — ~i''^'^ki ^0 = central spin 

k = 1, . . . , g, and g = n — 1. Note that, since the spec- 
tral polynomial Q2n{u) in equation H3.1(j|l is positively 
defined, the matrix lo is pure imaginary and, therefore, 
all frequencies fife are real (see Fig. 1). 

IV. DEGENERATE SOLUTIONS 

Here we identify a class of "resonant" nonstationary 
solutions for the BCS and central spin models with n 
degrees of freedom (n spins or pseudospins) but only m < 
n independent frequencies. 

While a typical dynamics of the BCS and central spin 
models with n (pseudo)spins is quasi-periodic with n 



independent frequencies, it is clear already on physical 
grounds that there must be solutions with fewer frequen- 
cies. For example, the only "spatial" symmetry of the 
BCS ground state is the U{1) symmetry of the order pa- 
rameter that translates into rotations of all pseudospins 
in Hamiltonians H2.1|l and (|2.t)|) around the z-axis. The 
corresponding trajectory in the phase space is a circle 
and, therefore, the motion is periodic with a single fre- 
quency (see also the discussion following equation H2.6() 
and below in this section). 

Generally, one can show'^^ that there are (resonant) 
tori with an arbitrary number to < n of independent 
frequencies in a vicinity of any point in the phase space. 
It should be noted though that the set of points for which 
m < n has a zero measure in the phase space just like 
the set of rational numbers on the real axis. Moreover, 
such points are typically unstable^ and provide seeds 
of chaotic behavior in integrable systems. This occurs 
because small perturbations are able to destroy resonant 
tori and let the corresponding trajectories escape to other 
regions of the energy shell. In contrast, the majority of 
invariant tori with n independent frequencies are only 
slightly deformed by perturbations. 

Interestingly, for the resonant (degenerate) cases con- 
sidered here we were able to reduce the BCS and central 
spin models with n spins to same models with only m 
spins. For example, solutions with to = 2 frequencies 
can be obtained from the BCS and central spin models 
for only two spins etc. Further, we will see that one spin 
solutions (periodic trajectories) are special among degen- 
erate solutions in that they correspond to the BCS energy 
spectrum. The solution that corresponds to the ground 
state is further exceptional in that it is stable against 
conservative perturbations because it minimizes the to- 
tal energy. In this case, even though there is only one 
frequency, the trajectory cannot leave the resonant torus 
because the latter coincides with the energy shell. The 
same applies to the solution that minimizes the central 
spin Hamiltonian. 

We start by determining conditions under which the 
number of independent frequencies is reduced. Since the 
frequencies are fixed by initial values of integrals of mo- 
tion, degeneracies occur only for special values of the inte- 
grals. On the other hand, a complete information about 
integrals of motion is encoded in the spectral polynomial 
Q2niu) defined in equations H3.3|l and (|3.4() . Specifically, 
we have seen in the previous section that the number of 
frequencies Vlk is equal to the number of branch cuts of 
the hyperelliptic curve z{u) = \J Q2n{u)- Evidently, the 
latter decreases by one when two roots of the polyno- 
mial Q2niu) coincide. Thus, as detailed below, merging 
2{n — m) roots of Q2n{u) we obtain solutions with m < n 
frequencies. 

Indeed, consider a situation when 2(n — to) roots co- 
incide, i.e. the spectral polynomial Q2n{u) in equation 
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(|3.3|) has the following special form: 

Q2n{u) = Q2,n{u) Jl {u - Ekf = Q2m (") W^^_,„ (w) 
k—m 

(4.1) 

We assume that double roots in equation H4.1|l are com- 
plex conjugate to each other so that coefficients of poly- 
nomials Q2m{u) and Wn-m{u) are real. Note that, since 
the polynomial Q2n{u) is positively defined (see the re- 
mark following equation H3.5|l ). so is the polynomial 
Q2m (u) ■ We will see below that degenerate solutions are 
completely determined by the "residual" spectral poly- 
nomial (52m (w). 

Now let us choose m — n separation variables Uj to 
coincide with double roots of the polynomial Q2n{u), i.e. 

Uj — Ej j — m, . . . ,n — I 

n-l n-l (4.2) 

W„_,„(u) ^ Y[{u- Ek) = Yl{u-u,) 

k=m j=m 

We observe that equations of motion for these variables 
are automatically satisfied because both sides of equa- 
tions H3.9|l and H3.10|l vanish. 

Next, we use equations (|4.1|l and H3.3|) to express eigen- 
values of matrix C in terms of the residual spectral poly- 
nomial Q2m{u) as follows 




where we expanded the ratio of polynomials Wn-miu) 
and Pn{u) into elementary fractions 

Wn-m{u) ^ Cg 

Pn{u) ' ' 

Since the degrees of polynomials Wn-m{u) and Pn{u) 
differ by m, there are m constraints on the values of C^, 
which can be derived by expanding equation H4.4|) in l/u. 
We have 

n-l 

^C,eJ-i=4m fc = l,...,m (4.5) 

9=0 

Values of Cq can be expressed through coefficients of 
Q2miu) by comparing residues at double poles at u = 
in equations H3.3|l and (|4.3|l . 

Cq = J^'' where Cq = ±1 (4.6) 

a\J Q2mi<^q) 

Note that equations l|4.5|l provide the following m con- 
straints 

i4-^L^ = = ctSkm k^l,...,m (4.7) 

9=0 J Q2m{£q) 



on 2m coefficients of a positively defined polynomial 
Q2m{u)- 

Values of integrals of motion Hq for which degenera- 
cies occur can be determined in terms of coefficients of 
the residual spectral polynomial Q2m(u) by comparing 
residues at simple poles at u — eq in equations (|3.3|l and 

631). 

Hq = E ' ^^l^H^ + ^C^qQM (4.8) 

In particular, the value of Hq - twice the energj*^ of the 
central spin model on degenerate solutions is 

n-l „ ~ 

fo = /iV^eoE^ + :^ (4.9) 

9=1 

where Sfe are coefficients of Q2m{u) defined as 

2m 

Q2rniu) = ^aku'' a2m = 1 (4.10) 

fe=0 

The BCS energy £bcs oc J2q ^q^q + const and the value 
of Jz (X J2q Hq can be obtained from the knowledge of 
Hq. However, a more convenient derivation is to compare 
expansions of equations ^6. '61 and H4.3|l in powers of l/u. 
We have 

n-l ~ 
Jz + a'^Cqe'^ = (4.11) 

9=0 
m — 1 

gSscs = - X! ^9(2e"'''^ + a2m-ie")- 

(4.12) 

— 2 

^2m-l , ~ 

— V aim-2 

So far, we determined initial conditions for which the 
number of independent frequencies for a system with n 
(pseudo)spins drops from n to m. We saw that these 
conditions are fixed by the residual spectral polynomial 
Q2m{u). Note that separation variables uj with j > m 
are also fixed by this polynomial because they are zeroes 
of equation H4.4|) . 

Thus, we are left with m dynamical variables - Uj with 
j = 1 , . . . , TO — 1 and J_ . Equations of motion (|3.9|) and 
(|3.1l)|l for these uj can be cast into a standard form (|3.11|) 
where n is now replaced with to and the original spectral 
polynomial Q2n{u) is replaced with the residual polyno- 
mial Q2m{u). On the other hand, it is clear that we 
would obtain similar equations of motion for to spins Sj 
governed by the BCS or central spin Hamiltonians. This 
analogy can be made precise if we identify the residual 
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polynomial Q2m{u) with the spectral polynomial of a sys- 
tem of m spins (cf. equation 1)3. i.e. 




3=0 



(4.13) 

where we replaced parameters tj with new ones 

m — 1 

Pm{u)=\{{u-l^) (4.14) 

g=0 

and chose eo = 0, consistent with the choice eo = (see 
below equation H3.10|l ). The rest of are arbitrary real 
numbers. 

BCS (|2.5|l and central spin (|2.9|) Hamiltonians for new 
classical spins Sj are 



m— 1 

Hbcs — ^ 2e 



m — 1 



Ho=J2 y^o • Sj - Bsl 7j 



g=0 

2 

eo - e J 



(4.15) 



Now, comparing equations of motion (|3.9|l and (|3.1U|) 
for m separation variables Uj for new spins with those 
for the remaining variables Uj for j = 1, . . . , m — 1 in the 
original problem with n spins, we see that they coincide if 
we re-scale the time variable for the central spin problem 
as follows 



t = _ t central spm 



(4.16) 



t = t 



BCS 



No re-scaling is necessary for the BCS model. Further, 
analyzing equations of motion for x and y components 
of the total spin (|3.9|) and 1)3. 10|) , we conclude that they 
coincide for original and new spins, i.e. 

J- = J- 

To express original n (pseudo)spins Sj{t) — 2(Kj(t)) in 
terms of m new spins Sj(f), we use equations 1)4.4(1 and 
1)4.2(1 to rewrite the matrix element B{u) given by equa- 
tion 1(3.7)1 as 



m — l 



n-1 



B{u) - J- n - E 



3 = 1 9=0 
n-1 



(4.17) 



B{u)P,n{u)Y, 



9=0 



u- eg 



Evaluating residues at u = €j in equation 14.17(1 . we ob- 
tain according to equation 1)3.8)1 



9=0 ^« 



(4.18) 



where constants Cj are given by equation ()4.6)) . Sim- 
ilarly, one can relate z-components of original spins to 
new ones 



sUt) = QP™(e,) 



m — l 

E: 

9=0 



(4.19) 



The solution for new spins (t) as functions of time can 
be read off directly from the general solution - equations 
l(X^ . and of the previous section. We 

only need to replace n with m, re-scale the time variable 
for the central spin model according to equation 1)4. 16)1 . 
and replace the spectral polynomial Q2n{u) with its re- 
duced counterpart Q2m{u) in the definition of matrices 
of periods ()3.16)l of hyperelliptic functions. 

For future references, let us also write down equations 
for components of original spins Sj{t) — 2{'Kj{t)) and 
evolution of J-{t) for degenerate solutions in terms of 
m variables Uj = Uj. Evaluating residues at w = e^ in 
equation 1)4.17)1 and using equation ()3.8)l . we find 



(4.20) 



With the help of equation 1)4.1)) and 1)4. 20)1 . evolution 
equations ()3.9)) and 1)3.10)1 for J_ become 



din J_ 



(-l)"ieo 



/ao 



Uj central spin 



(4.21) 



a2m-l + 2j2uA BCS 



We see that degenerate solutions with m < n frequencies 
can be parameterized by m auxiliary spins both for the 
BCS and central spin model. Therefore, corresponding 
trajectories live on an m-dimensional torus. By sym- 
metry, for the same initial conditions, trajectories of all 
Hamiltonians Hq in system ()2.8)l live on the same torus. 
This implies, for example, that a slight change of ini- 
tial values of integrals of motion 14.8)1 will induce small 
oscillations along the remaining n ~ m directions. 

Let us consider several examples of degenerate solu- 
tions. 

One spin solutions, m = 1. These solutions reproduce 
the BCS solution for the ground state and excitations of 
the BCS model (see also the discussion following equa- 
tion 

Since the residual spectral polynomial is positively de- 
fined, we can parameterize it as Q2(u) = {u — /i)^ + Aq. 
For m = 1, there is only one condition in equation 1)4.7)1 



n-1 



(4.22) 



For the BCS model a = 2/g and we recognize equa- 
tion ()4.22)) as the BCS gap equation. Similarly, equation 
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(|4.11() becomes the BCS equation for the chemical poten- 
tial. 

The energy of BCS model H4.12|l for m = 1 is 



n-l 



3=0 



We see that a choice of signs Cj — 1 for all j yields 
the BCS ground state energy, while if we choose one 
of Ej to be negative, eq — —1 and Cj = 1 for j ^ g, 
the solution corresponds to a pair excitation of energy 
2yJ {cq — + Aq. Similarly, one obtains solutions with 
two and more excitations by choosing several signs Cj to 
be negative. 

According to equations (|4.15|l . the BCS and central 
spin models for m = 1 reduce to the following single spin 
models: 



Hq — —Bsq Hbcs 



(4.24) 



Then, equation l|4.13|l provides an alternative 
parametrization of the residual spectral polynomial 

Q2{u) 



(4.25) 



The solution of equations of motion is 

J_ (t) = (i) = aAoe-"3*+''^ (4.26) 

where the frequency /3 is different for BCS {(3bcs) and 
central spin (/3o) models 



Pbcs = S^i (3o 



eo = ±l 



Components of original spins s.j{t) — 2(Kj(f)) are given 
by equations H4.18|l and (|4.19(l . We have 



(4.27) 

Note that one spin solutions constructed above do not 
minimize central spin Hamiltonian H2.9|l . Instead, we ob- 
serve that in this case the configuration of spins with 
minimum energy is s| = — sgn7j. Similarly, we see from 
equations H2.8|l that configurations that minimize Hamil- 
tonians Hq are s| = — sgn(eg — fj). Interestingly, accord- 
ing to the definition of pseudospins l|2.2|) . one of these 
configurations is the unperturbed Fermi ground state. 
Since all spins in these configurations are oriented along 
the z-axis, they are also stationary with respect to BCS 
Hamiltonian H2.5|l . 

Two spin solutions, m = 2. In this case there are 
two constraints H4.7|) on four coefficients of a positively 



defined polynomial Q^^u). 



E 



E 



q=0 jQ4eq) q=0 \ Qi{eq) 



a Cq — ±1 



4 

E 

fc=0 



aku 



a4 = 1 



Equations H4.9|l and (|4.12|l show that the energy of these 
solutions can take arbitrary values in the range of ener- 
gies allowed for the BCS and central spin models. 
The equivalent two spin problems (|4.15() are 

i^BCS = 2?isJ - I J+J_ ^ llso -Si - (4.28) 

We are left with only one separation variable ui = ui 
and equations H3.11|l now read 



ul+AQ4ui)^0 BCS 

otQuI + Q4{ui) = central spin 



(4.29) 



while for the total (pseudo)spin, using equation 14.21|l . 
we obtain 



din J- i . 

— — ±—^=ui central spm 

at V ao 



(4.30) 



dlnJ_ 



ias + 2iui BCS 



Components of individual spins can be expressed either 
through ui and J_ or in terms of new spins using equa- 
tions and (j^TH|) 

sj = CjJ-iej - wi) = Cj(ej -?i)so (?) + CjejS^(t) 

(4.31) 

A special case of two spin solutions when the four coef- 
ficients of a positively defined polynomial Qi^u) can be 
parameterized by three numbers, Q^iu) = {{u — lu)"^ + 
-I- A^)^ — 4A^ A^, has been previously discovered in 
Ref. 0. The general two spin solution is obtained from 
equations H4.29I4.31() . which are solved by setting n — 2 
(g = 1) in equations (|TT7|) . and where hy- 

perelliptic functions now become ordinary elliptic func- 
tions. 

Similarly, one can obtain three, four etc. spin solutions 
(degenerate solutions with m = 3, 4 etc. frequencies) in 
terms of hyperelliptic functions of genus g = 2,3 etc. 
However, as discussed in the beginning of this section, 
solutions with m < n frequencies are nonrepresentative 
of the dynamics anywhere in the phase space and are 
expected to be unstable in the sense of KAM theorem. 

Moreover, we emphasize that degenerate solutions de- 
rived here are not all solutions with m < n indepen- 
dent frequencies. Indeed, it is clear that for a generic 
point in the phase space all roots of the spectral polyno- 
mial Q2n{u) are distinct and, by continuity, the distance 
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to the closest point where two roots coincide is finite. 
On the other hand, there are points with any number 
of frequencies arbitrarily close to any point in the phase 
space. Other solutions with m < n frequencies can be 
constructed using the reduction theory for hyperelliptic 
functions^. 

Note also that degenerate solutions with m frequencies 
contain ones with fewer frequencies and can be further 
reduced to m — 1, to — 2 etc. Geometrically, they live on 
a 2TO-dimensional surface in the 2n-dimensional phase 
space (see the remark below equation (|4.21|l '). Surfaces 
with smaller m are embedded into ones with larger to. 
Interestingly, periodic trajectories that correspond to the 
BCS energy spectrum (m — 1) belong to all these sur- 
faces. 



first a simpler setup of a single dot connected to two po- 
larized leads. This setup leada^ to an integrable model 
very similar to the central spin model. 

Further, it might be possible to use the BCS model 
(|2.t)|) with few classical spins to describe a number of 
grains connected to each other by Josephson junctions. 

Other interesting applications include pairing phenom- 
ena in cold fermion gasesi: modeled either by BCS or 
Dicke model and in superconducting metals (see e.g. the 
discussion in Ref . Issj) . In these cases, a careful identi- 
fication of observable robust features of the solution is 
needed. 

An interesting open problem is the evaluation of lead- 
ing finite size (quantum) corrections to the general solu- 
tion obtained in this paper. 



V. APPLICATIONS AND OPEN PROBLEMS 

Our results can be directly applied to the problem of 
decoherence within the central spin model. In this case, 
we believe, a complete exact answer can be obtained us- 
ing asymptotics of hyperelliptic functions. 

Another possible application is to experiments^ on 
electron transport in spin blockaded semiconductor dou- 
ble dots. In this connection, it might be useful to analyze 
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